data <- readRDS('data_sim_p_1000.rds')
for(x in names(data)) assign(x, data[[x]], envir = globalenv())
U_full <- data$U
rm(data)
S.data.time <- S.data$obs
S.data.time.log.sum <- sum(log(S.data.time))
model <- SAEM_model(
function(sigma2, ...) -n/(2*sigma2),
function(phi1, phi2, phi3, ...) mean((Y - get_obs(phi1, phi2, phi3) )^2 ), 'sigma2',
# === Variable Latente === #
latent_vars = list(
# === Non linear model === #
latent_variable('phi', dim = G, size = 3, prior = list(mean = 'mu', variance = 'omega2'),
add_on = c('zeta(phi1 = phi1, phi2 = phi2, phi3 = phi3, ...)' )),
# === S.data model === #
latent_variable('b', prior = list(mean = 'barb', variance.hyper = 'sigma2_b'),
add_on = c('zeta(b = b, ...) +',
'sum(h$eval(b = b, ..., i = c(1,2)))' )),
latent_variable('alpha', prior = list(mean = 'baralpha', variance.hyper = 'sigma2_alpha'),
add_on = c('zeta(alpha = alpha, ...) +',
'alpha*h$eval(alpha = alpha,..., i = 3)'))
),
# === Paramètre de regression === #
regression.parameter = list(
regression_parameter('beta', 1, function(...) SPGD(5, theta0 = beta,
step = 0.05, lambda = 1/sqrt(N),
normalized.grad = T,
zeta.der.B, N, zeta.B,
Z$alpha, Z$phi1, Z$phi2, Z$phi3,Z$b) )
)
)
# --- Initialisation des paramètres --- #
parameter0 <- parameter %>% sapply(function(x) x* runif(1, 1.1,1.4))
parameter0$barb <- 10
parameter0$sigma2 <- 1
parameter0$baralpha <- 1
#===============================================#
load.SAEM(model)
oracle <- maximisation(1, do.call(S$eval, var.true), parameter, var.true)
#==============================================================================#
init.options <- list(x0 = list(phi = c(1,80,4), b = parameter0$barb, alpha = parameter0$baralpha),
sd = list(phi = c(.05, 1.5, .5), b = 1, alpha = .1) )
SAEM.options <- list(niter = 200, sim.iter = 5, burnin = 190,
adptative.sd = 0.6)
saem
p <- 4
U <- U_full[,1:p] #Mise à jours des cov
parameter$beta <- parameter$beta[1:p] #Mise à jours de beta
parameter0$beta <- runif(p, min = -1, max = 1) #initialisation de beta
SAEM.options$RDS <- paste0(params$rds_filename, '_', p)
res <- run(model, parameter0, init.options, SAEM.options, verbatim = 3)
## [1] "SAEM execution time = 45min 44sec"
## $plot_parameter

##
## $plot_MCMC

##
## $plot_acceptation

## [[1]]

##
## [[2]]

##
## [[3]]
Result of the SAEM-MCMC
|
|
sigma2
|
mu.1
|
mu.2
|
mu.3
|
omega2.1
|
omega2.2
|
omega2.3
|
barb
|
baralpha
|
|
Real value
|
0.0025
|
0.9032
|
89.9575
|
5.0079
|
0.0039
|
35.9179
|
0.6948
|
30.0000
|
45.0000
|
|
Estimated value
|
0.0025
|
0.8979
|
90.0274
|
4.9032
|
0.0040
|
34.2305
|
0.5823
|
18.2806
|
30.2388
|
|
Rrmse
|
0.0101
|
0.0059
|
0.0008
|
0.0209
|
0.0098
|
0.0470
|
0.1619
|
0.3906
|
0.3280
|
p <- 20
U <- U_full[,1:p] #Mise à jours des cov
parameter$beta <- parameter$beta[1:p] #Mise à jours de beta
parameter0$beta <- runif(p, min = -1, max = 1) #initialisation de beta
SAEM.options$RDS <- paste0(params$rds_filename, '_', p)
res <- run(model, parameter0, init.options, SAEM.options, verbatim = 3)
## [1] "SAEM execution time = 46min 29sec"
## $plot_parameter

##
## $plot_MCMC

##
## $plot_acceptation

## [[1]]

##
## [[2]]

##
## [[3]]
Result of the SAEM-MCMC
|
|
sigma2
|
mu.1
|
mu.2
|
mu.3
|
omega2.1
|
omega2.2
|
omega2.3
|
barb
|
baralpha
|
|
Real value
|
0.0025
|
0.9032
|
89.9575
|
5.0079
|
0.0039
|
35.9179
|
0.6948
|
30.0000
|
45.0000
|
|
Estimated value
|
0.0025
|
0.9004
|
90.1263
|
4.9105
|
0.0044
|
34.6121
|
0.5214
|
17.7837
|
30.0382
|
|
Rrmse
|
0.0136
|
0.0031
|
0.0019
|
0.0194
|
0.1135
|
0.0364
|
0.2495
|
0.4072
|
0.3325
|
p <- 100
U <- U_full[,1:p] #Mise à jours des cov
parameter$beta <- parameter$beta[1:p] #Mise à jours de beta
parameter0$beta <- runif(p, min = -1, max = 1) #initialisation de beta
SAEM.options$RDS <- paste0(params$rds_filename, '_', p)
res <- run(model, parameter0, init.options, SAEM.options, verbatim = 3)
## [1] "SAEM execution time = 46min 12sec"
## $plot_parameter

##
## $plot_MCMC

##
## $plot_acceptation

## [[1]]

##
## [[2]]

##
## [[3]]
Result of the SAEM-MCMC
|
|
sigma2
|
mu.1
|
mu.2
|
mu.3
|
omega2.1
|
omega2.2
|
omega2.3
|
barb
|
baralpha
|
|
Real value
|
0.0025
|
0.9032
|
89.9575
|
5.0079
|
0.0039
|
35.9179
|
0.6948
|
30.0000
|
45.0000
|
|
Estimated value
|
0.0026
|
0.9025
|
90.1870
|
4.9241
|
0.0047
|
34.7273
|
0.5773
|
15.3291
|
27.0420
|
|
Rrmse
|
0.0189
|
0.0007
|
0.0026
|
0.0167
|
0.1872
|
0.0331
|
0.1690
|
0.4890
|
0.3991
|
p <- 1000
U <- U_full[,1:p] #Mise à jours des cov
parameter$beta <- parameter$beta[1:p] #Mise à jours de beta
parameter0$beta <- runif(p, min = -1, max = 1) #initialisation de beta
SAEM.options$RDS <- paste0(params$rds_filename, '_', p)
res <- run(model, parameter0, init.options, SAEM.options, verbatim = 3)
## [1] "SAEM execution time = 36min 26sec"
## $plot_parameter

##
## $plot_MCMC

##
## $plot_acceptation

## [[1]]

##
## [[2]]

##
## [[3]]
Result of the SAEM-MCMC
|
|
sigma2
|
mu.1
|
mu.2
|
mu.3
|
omega2.1
|
omega2.2
|
omega2.3
|
barb
|
baralpha
|
|
Real value
|
0.0025
|
0.9032
|
89.9575
|
5.0079
|
0.0039
|
35.9179
|
0.6948
|
30.0000
|
45.0000
|
|
Estimated value
|
0.0025
|
0.9037
|
90.0198
|
5.0017
|
0.0040
|
35.3809
|
0.6134
|
4.2174
|
8.2451
|
|
Rrmse
|
0.0046
|
0.0005
|
0.0007
|
0.0012
|
0.0277
|
0.0150
|
0.1172
|
0.8594
|
0.8168
|